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ABSTRACT 



The problem of obtaining a control policy which is optimum in some 
respect for a specified system is one of the more pressing problems of 
automatic control theory today. L. S. Pontryagin stated his Maximum 
Principle in 1956 and L. I. Rozonoer discussed and extended the work in 
a series of articles published in 1959. 

The Maximum Principle, as it applies to the ’’free right end" problem 
for a non-linear, one-dimensional system, is utilized in conjunction with 
the CDC 1604 digital computer to obtain optimal control policies for 
several different cost functions. 

The authors wish to express their appreciation for the encourage- 
ment given them by Dr. Harold A. Titus of the U„ S. Naval Postgraduate 
School and to express their gratitude to Miss Mary E, Haynes of the 
Computer Facility, U. S. Naval Postgraduate School whose patient program- 
ming assistance materially aided in the successful completion of this 
thesis. 
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1 . Introduction, 



One of the most important problems of automatic control theory is 
the problem of creating systems which are optimal in some prescribed 
sense. Since the advent of satellites and space travel, a set of problems 
of the type to be discussed here have taken on greater importance. Since 
many of the most interesting problems in this field have system equations 
for which the methods of the Calculus of Variations fail, new mathemati- 
cal methods are required for their solution. 

Suppose that the system to be investigated may be described by an 
equation of the form 



ff = G(X,U) . X(0) = C (1-1) 

where X(t) is the state variable and U(t) is the control variable.* The 
problem then is to determine U(t) such that the system is optimized 
according to some specified criterion or cost function. 

Two types of process are of particular importance: 

(1) A desire for X(t) to remain as close to a prescribed path as 
possible throughout the duration of the process and, 

(2) A desire for the final value of X(t) to be a prescribed value 
or in a prescribed state (Terminal Control). 



v Eqn. (1-1) is written in one dimension and this investigation will 
be confined to systems of the first order. However, all methods and 
techniques employed may easily be extended to n-th order systems where 
X and U would be replaced by the vectors X and U. 
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The type of problem as set forth in case (1) is the type investigated 
here and is the ’’Free Right End" problem as described by L. I. Rozonoer 
in Section I of his paper dealing with Pontryagin's Maximum Principle / 1/ . 

A non-linear system was hypothesised and methods of obtaining a 
control policy which would provide an optimum system according to some 
specified criterion were investigated. 

The system chosen for the investigation may be described by a non- 
linear differential equation of the first order which has the following 
functional form: 

X = f(x,u) , X(0) = C (1-2) 

in which one may consider X as representing the error or deviation from 
a prescribed path and U as representing the control effort applied. The 
problem was to determine the control policy U, such that some criterion 
function of the form 

J = ^G(X,u)dt (1-3) 

would be minimized. Such functions as described by Eqn. (1-3) with the 
constraint of Eqn. (1-2), which evaluate the performance of a system are 
commonly referred to in the literature as "criterion" or "cost" functions. 
Sometimes they can be expressed simply in terms of integrals. At best, 
the choice of a cost function J, is a compromise between a desired crite- 
rion of goodness for the control design and one which leads to a more 
tractable mathematical analysis. 

The specific system investigated is described by the first order, 
non-linear differential equation 

X = AX + B Sin + CU (1-4) 
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with 


the constants having the 


f o 1 lowing values^ 






A = -0.1 


B « +1.0 


C » +.25 




Several different cost functions, subject to the 


constraint of 


Eqn. 


(1-4) were proposed for 


the investigations 






J = 


C (X 2 +U 2 )dt 
' 0 


(1-5) 




J = 


( (x 2 )dt 
' 0 


(1-6) 




J = 


f (u 2 )dt 


(1-7) 



The problem of obtaining a control policy to effect the minimization 
of these criteria (each of which leads to a different set of adjoint 
equations) will be treated in detail below. 

Qualitatively, minimizing Eqn. (1-5) would correspond to minimizing 
the error and the control effort to the system. Minimizing Eqn. (1-6) 
would correspond to minimizing only the error, while minimizing Eqn. (1-7) 
would minimize only the control effort of the system with no regard for 
the error. 

As mentioned above, this investigation deals with the "Free Right 
End" problem. As a result, the time interval over which the integrals, 
Eqs. (1-5), (1-6) and (1-7), were evaluated was fixed at 10 seconds. 

By the very nature of the formulation of the problem (minimization 
of the integral J = j G(X,u)dt with the constraint X * f (X,U), having 

X(0) specified, but not specifying the value of X at the final time T), it 
falls within the province of the Calculus of Variations. Occasionally 
such classical methods can be employed to determine the optimal control 
policies. 

The authors however, attempted to apply the work of a Russian control 
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theorist, L. S. Pontryagin, as it is presented in a series of articles 
by L. I. Rozonoer /1,2,3/. With the aid of these papers and the CDC 
1604 digital computer, satisfactory control policies were obtained. 

Briefly recapitulating, the purpose of this thesis is to use the 
digital computer, employing various programming methods, to search for , 
compute, and design an optimal control which will satisfy the various 
criteria (cost functions) as applied to Eqn. (1-4 ) . 
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2. The uncontrolled system . 

The uncontrolled system as described by the following equation 

X = -0.1X + Sin X C 2 ' 1 ) 

was investigated first. A knowledge of the behavior of the system with- 
out restraint or control applied was deemed necessary in order that 
investigations of the system with control could be properly interpreted 
and understood. 

The Donner Analog computer (Model 3100) was employed in order to 
obtain the phase plane plot of velocity (X) versus trajectory (X). Fig. 
(2-1) shows the computer diagram used to obtain the phase plane desired. 







FIG. (2-1). Computer diagram utilized in order to obtain the 
Phase Plane (X vs.X) for the uncontrolled system. 

As no sine function generator was available, the sine function was 

simulated as shown within the dotted rectangle.* 



*The sine function was simulated by constructing the analog of the 
second order differential equation X + X = 0 whose solution is known to 
be X * sin(t). 
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The output 



-Y = 20 Sin t -2t (2-2) 

was recorded on an X-Y Plotter utilizing the real time variable input to 
the X coordinate. Fig. (2-2) shows the results of this investigations, 
the phase plane plot for the uncontrolled system. 




FIG. (2-2). Phase Plane Plot of the uncontrolled system showing 
Velocity (X) on the ordinate plotted against Trajectory (X) on 
the abscissa. 

Since the computer output was in terms of Y and Time, it was neces- 
sary to convert these coordinates to the desired coordinates X and X. 

In order to scale the abscissa it was necessary to find the point on 
the curve corresponding to n radians. Since the curve as described by 
Eqn. (2-2) is a linear combination of a sine curve and a straight line 

U = -Kt (2-3) 
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the points where the resultant curve intersects the straight line are 
multiples of it. The points on the abscissa lying directly above these 
points of intersection are, therefore, also multiples of it radians. In 
this manner the X coordinate was transformed from a time scale to a 
distance scale. (Since only qualitative results as to velocity were 
desired the Y *xis was not rescaled and no transformation of coordinates 
was made). 

Investigation of Fig. (2-2) yielded the following interesting informa 
tion concerning the existance of equilibrium points: There are two stable 

equilibrium points located at X = 2.84+ and at X - 8.41+ and there are two 
unstable equilibrium points located at X = 0.0 and at X = 7.02+. Thus, 
the uncontrolled system may be expected to move toward one of the stable 
equilibrium points depending on the initial value of X as indicated by 
the arrows on the curve of Fig. (2-2). 

In addition to the above investigation, the GDC 1604 digital computer 
was employed as an additional method of determining the equilibrium points 
By choosing selected values of initial X (both positive and negative) and 
solving the differential equation (by the Runge-Kutta method), equilibrium 
points were obtained which agreed favorably with those obtained by the 
graphical analysis of the analog solution. (Fig. (1-2) and Fig. (1-3) 
in Appendix I are the computer graph solutions obtained from this invest i= 
gation. Tabular output also was obtained but is not included here due 
to the volume obtained and due to the relative unimportance of this sort 
of detail). 

In order to obtain values of the equilibrium points to an accuracy 
greater than that possible by either of the preceeding two methods, 
standard mathematical methods were utilized. The derivative was set 
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equal to zero and the values of the equilibrium points were obtained to 
three decimal point accuracy with the aid of sine tables, a desk cal- 
culator and iterative proceedures. The results obtained again corres- 
ponded to and verified the earlier results, yielding stable equilibrium 
points at X = 2.852+ and at X = 8.416+ with two unstable equilibrium 
points at X = 0.00 and X = 7.068+. Only positive equilibrium points 
were calculated as this investigation was to be limited to positive 
values of X. 

Interpreting the significance of the sum of these investigations, 
it is possible to define the regions of stability and instability for the 
uncontrolled system as shown below in Fig. (2-3) 



00 

8 

6 

4 

2 



Stable equilibrium point 



8.416+ 



1 

7.068+ 



2.852+ 



Un stable e q ui li brium p oint 



Stable equilibrium point 



0 Un stable eq ui 



< 



n 

55 



Stable equilibrium p oint 



-4i 



Pnst a ble equilib rJL u rn point 



Stable equilibrium poin t 



-2.85+ 



-7.07+ 

-8.52+ 



00 

FIG. (2-3) Regions of stability and instability for the 
uncontrolled system. 
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3. Maximum principle in optimal system theory . 



As stated above, in an effort to obtain a solution to the problem 
of obtaining optimal controls, the authors utilized the "Maximum 
Principle" as hypothesized in 1956 by L. S. Pontryagin on the basis of 
work performed by him, V. G. Boltyanskii and R. V. Gamkrelidze. / 4/ 
Pontryagin' s Maximum Principle is presented in some detail by L. I. 
Rozonoer in a series of articles appearing in "Automatika i Telemakhanika" 
in 1959. ("Automation and Remote Control" presents an English transla- 
tion of this Russian Journal. /1.2.3/). A brief resume of the most 
important concepts will be made here in order that a common ground for 
further discussion may be establ shed. 

Rozonoer' s papers deal with the problem of obtaining a control which 
is optimum in some sense for a system which may be described by a set of 
differential equations of the n-th order:* 



where X^, ..., X^ are the parameters of the system and lL,...,U r are 
the positions of the controlling elements. 

Rozonoer shows in his paper that the problem of optimizing the 
system with respect to an integral 



leads to the problem of optimizing with respect to coordinates.** At this 



*Any n-th order differential equation may be expressed in terms of 
n first order differential equations involving n variables. 

**L. I. Rozonoer, L. S. Pontryagin* s Maximum Principle in the Theory 
of Optimum Systems, "Automation and Remote Control", Vol 20, pl291. 



X A = f(X i( U i( t) i=l,...,n 



(3-1) 




(3-2) 
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point an additional variable may be introduced 



X 



n+1 




which allows one more differential relation 




(3-4) 



to be added to the system specified in Eqn. (3-1). The problem thus 
becomes one of optimizing the n+lst system coordinate at the final 
moment of time. 

Specifically, the problem of optimizing a linear function of the 
final values of all the coordinates of the system, that is, the quantity 



where are certain constants, is developed in detail in Rozonoer's 
paper. In the discussion of this problem, the theory is first developed 
for the case in which the right end of the trajectory, X(t), is not fixed. 
That is, there are no restrictions imposed on the final values of the 
coordinates .* 

A variable vector, P(t) = P^(t),..., P^(t), which has a direction 
at time t = T opposite to the direction of the vector C = is 

introduced at this point 



It is assumed that the moduli! of the vectors P(T) and C are equal. 



*For the entire problem considered here, the final coordinate was 
left free. Only the duration of the problem was fixed at 10 seconds. 




k=l 



(3-5) 



P i (T) = -°i 1 = 1 



(3-6) 
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The variables P.(t) are subject to a set of differential equations 



n k ;t) 

,(t) = - ^TP o -i-i iLJL £ ,• ; i = 

1 



3x. 

1 



(3-7) 



(Let us note here that if any control, U(t), is given, the vector P(t) 
is uniquely determined from Eqs. (3-7), where conditions Eqs. (3-6) 
play the role of boundary conditions.) From Eqs. (3-5) and (3-6) bound- 
ary conditions for the final values of the adjoint variables, P , may be 
obtained. 

Rozonoer continues, and points out that if the following expression 
is formed. 



h(x 



,p,u,t) = p A 



(3-8) 



S=1 



that Eqs. (3-2) and (3-7) may be written in the form 



X 

\ 



p. 

1 



c)H 

-Sx 4 



i = i,. 



• > n 



(3-9) 



Since the X^(0) are specified in the problem statement and the 
P^(T) may be found from Eqs. (3-6) and (3-5), the boundary conditions 
on Eqn . (3-9) are 



X i (0) = X° , P^T) = - C ± i = (3-10) 

The H-function is analogous to the Hamiltonian in analytical 
mechanics, and the vector P(t) to the impulse vector. Rozonoer proves 
in his paper that, according to the maximum principle, the H-function 
must be maximum (minimum) in U for any values of X and P in order to 
obtain the optimum condition. 

Now well known principles of Calculus are resorted to and the 
partial derivative of the H-function is taken with respect to U and set 
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equal to zero. 



= 0 (3-H) 

Eqn. (3-11) yields a relation between the control U and the adjoint 
variable P. 

Since the cost function to be minimized has been defined as a 
linear combination of the final values of the coordinates of the object 
system, Eqn. (3-5), and this function added to the system as the n+lst 
coordinate, it may be seen from Eqn. (3-9) that 




P , = 0 

n+1 



(3-12) 



in all cases. This allows immediate integration of P . yielding 

n+l 



P , = CONSTANT (3-13) 

n+i 

Since P^ + ^ Is constant and we know from Eqn. (3-10) that p r+ j00 » “ C n+ i 
we see immediately that 



" C n + 1 (3 ‘ 14) 

Also note that will always be equal to one and that 

C. = 0 i« 1,2, ...,n 

Thus we see that p n+ ^ = “1 which allows some simplification in Eqs. (3-9). 

As a result of the above manipulations, 2n-l differential equations 
are obtained, the solution of which yield the desired optimum control. 
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4. The cost function J = 




(x 2 +u 2 -) dt. 



The first system investigated was one in which the integral 



J = 



(X 2 +U 2 )dt 



(4-1) 



'0 



was to be minimized. The system to be controlled was described by the 
differential equation 



for which the constants A, B, and C had the same values as specified in 
Section 2; namely, -0.1, 1.0, and 0.25 respectively. The time interval 
T, was fixed at ten seconds. (Some investigation as to the effect of 
extending the time interval to 20 seconds was investigated and is dis- 
cussed below.) 

Since the investigation of the uncontrolled system showed several 
points of stable and unstable equilibrium for initial values of X 
between zero and 10, an initial value of X equal to 20 was chosen in 
order to include the effects of the equilibrium points. 

The problem was approached according to the principles as set forth 
in Rozonoer's paper, / 1 / , dealing with Pontryagin’s Maximum Principle. 

In order to apply Pontryagin's Maximum Principle, it is necessary 
that the system be described by an equation of the form specified by 
Eqn. (3-1); i.e., a first order differential equation, either vector or 
scalar. Since the specific system under investigation, Eqn. (4-2), is of 
first order, no additional manipulation was required. 

Then, in accordance with Rozonoer's statement that optimizing with 
respect to an integral leads to the problem of optimizing with respect to 



X = AX + B SinX + CU 



(4-2) 
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the final values of the coordinates, the cost function was added as 
the n+lst coordinate 




The system equations now have the form 

X x = AX X + B SinX 1 + CU 

x 2 = X^ + u 2 



(4-3) 



(4-4) 



Next, forming the final value functional S, Eqn. (3-5), 
S 



= ^C.X.(T)= C 2 X 2 (T) = J^+U^dt 



(4-5) 



the values of the constants C. are obtained, 

i 9 




(4-6) 



and, since we know from Eqn. (3-6) that the adjoint variables are equal 
to but of opposite sense to the at the final time T, we now have the 
boundary conditions 



P X (T) = 0 

p 2 (t) = -1 



(4-7) 



Forming the H-Function, Eqn. (3-9), 



H(P,X,U,t) = ^Tp^ = 



(4-8) 



and making appropriate substitutions for the the H-function becomes 



H = + B SiriXj^ + CO) + P 2 (X 2 + U 2 ) . 



(4-9) 



From Eqs. (3-9) and (4-9) it is now possible to obtain the adjoint 
equations below: 
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P 1 = -(A + B Cos): i )P 1 - 2X 1 P 2 



(4-10) 



= 0 



Applying the principles of Calculus to the H-function in order to 
obtain the minimum with respect to the control U, a relation between 
the adjoint variables Pi and the control variables is found. 



Ait 

If = CP X + 2UP 2 = 0 



CP 1 
U ~ 2P_ 



(4-U) 



Noting that the equation for P 2 , Eqn. (4-10), is readily integrable, 
it is seen that P 2 (t) is equal to a constant. Since P 2 (T) = -1, Eqn. 
(4-7), it follows that P 2 must equal -1 at all times. Knowing the value 
of P 2 (t), it is possible to simplify Eqn. (4-11). 






(4-lla) 



Finally, making the appropriate substitutions for the constants 
A, B, and C and for the adjoint variable P 2 (t) and replacing U by P^/8, 
the following set of differential equations are obtained: 



X x = -.1^ + Sin + 


.03125 P x 


(4-12) 






(4-13) 


^ = -(.1 + Cos X 1 )P 1 


+ 2X X 


(4-14) 



The solution of the above set of equations is the problem remaining at 
this point. 

Numerical methods and the CDC 1604 digital computer were utilized 
in order to obtain a solution to this problem. Specifically, the Runge- 
Kutta method of evaluating first order differential equations, as program- 
med for the CDC 1604, was employed. (See Appendix V). 
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In order to obtain solutions to the set of differential equations 
by numerical methods, initial conditions are required. Since initial 
conditions are known only for the X variables, (X^(0)=20, X2(0)=0), and 
not for the adjoint variable ?£, the immediate problem becomes one of 
searching the adjoint space for the initial conditions which will yield 
the desired solution.* 

The first method proposed in an effort to obtain P*(0) was one in 
which several P(0) would be selected in an effort to bracket P*(0). Then, 
through some selective iterative scheme, the bracket size would be reduced 
and, eventually, the P*(0) yielding THE OPTIMUM system could be obtained. 
(This method of solution is commonly called a "hill climbing" technique). 
There was, however, no evidence that the variation of the cost function, 

J, versus the initial values of P would be smooth, and the possibility of 
"homing in" on some local minimum rather than the true minimum was present. 
(The investigation of the system showed that local minima did exist in 
many cases.) 

In view of the above, the authors decided on another approach to the 
problem of obtaining P*(0). A reasonable range of values for P*(0) was 
guessed. Then, the differential equations for the system were evaluated 
for X^(0) = 20, X2(0) = 0, and P^(0) ranging from -250 to +50. This 
investigation yielded values of cost function for 300 integral values of 
P(0) as shown in Fig. (4-1). 

A brief comment on the details of programming this problem might be 



*It should be pointed out here that any solution generated from Eqs. 
(4-12), (4-13) and (4-14) will be optimal with respect to the initial 
conditions chosen. However, of all these optimal solutions, there is one 
"best" solution. The P(0) which yields THE OPTIMUM solution will be 
designated P*(0). (The asterisk when applied to any other variable will 
likewise refer to THE OPTIMUM system.) 
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40000 




FIG. (4-1). Cost function vs. Initial value of P. 

X(0) =20 J = }(X 2 +U 2 )dt 

appropriate here. Since it had been decided to obtain a P*(0) by 

calculating the cost function for several P(0), the next question to 

resolve was in regard to the range of values of P(0) to investigate. It 

was suggested that a first guess might be that the initial control, U(0), 

should be at least as large as X(0) but of opposite sign. Utilizing Eqn. 

(4-lla) and recalling that X(0) = 20, a first guess puts P(0) at about -160. 

In an effort to allow for some uncertainty and also to obtain information 

as to the behavior of the cost function for P(0) not in the vicinity of 

P*(0), a range of F(0) was selected from -250 to +50 as stated above. 

Fig. (4-1) was the result of this initial investigation. (Not shown in 

Fig. (4-1), but obtained in a tabular output were the final values of the 

trajectory, X(T).) 

Locating the P*(0) which yields the minimum cost function and then 
repeating the above investigation in the neighborhood of the first deter- 
mined P*(0) +€ allows P*(0) to be determined to any degree of accuracy 
desired. 
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The major disadvantage of the above method is the time required to 
obtain P*(0) to a reasonable accuracy. In excess of 11 minutes was re- 
quired for the above investigation. (It was not possible to obtain tra- 
jectories or control policies during this first investigation due to 
storage limitations of the computer.) At that, the resulting P*(0) was 
accurate only to an integral value. 




FIG. (4-2). Optimal control policies, U(t), showing variation 
with the accuracy to which P*(0) is computed for J = J(X 2 +U 2 )dt. 

It should be noted that the accuracy to which P*(0) was determined 

affected both the control policy and the trajectory. Fig. (4-2) shows 

the effect on control that the accuracy of the P*(0) had. 

The OPTIMAL Control Policy obtained from the above computations 
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yielded a cost function of 1099.52 versus a cost function of 1978.09 for 
the uncontrolled system, an improvement* of 44.8%. 




FIG. (4-3). Trajectory vs. Time for the OPTIMAL system (solid 
curve) and the uncontrolled system (dotted curve). 

Fig. (4-3) shows the trajectories obtained for the uncontrolled 

system and the optimally controlled system. Note that the trajectory 

obtained for the optimal case as shown in Fig. (4-3) appears to approach 

a constant value of about 2.7 rather than a final value of zero. However, 

recall that the cost function to-be minimized here is a combination of 

error, X, and control effort, U, and also that the uncontrolled system 

had a stable equilibrium point at X = 2.84. One might assume, therefore, 



^Percent improvement is defined as 

Percent improvement = J (uncontrolled) - J (controlled ) 

J (uncontrolled) 
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CONTROL 



that the effort necessary to reduce the error below this value would 



exceed any reduction in the cost function that would result from such a 

TIME (Sec) 

reduction in X. 




FIG. (4-4). Optimal control policy for X(0) = 20 showing 
local maxima and minima. 

Fig. (4-4) shows the optimal control policy obtained for this system 
with X(0) = 20. It is very interesting to note that the points of local 
maxima and minima occur when the trajectory. Fig. (4-3), has a value such 
that the velocity, X, of the system is also very close to a maximum or 
minimum. (See Fig. (2-2) which is the velocity-displacement phase plane 
for the uncontrolled system). 

Having investigated the system with the initial value of X set at 
20, the next part of the investigation involved the determination of the 
effect of the initial value of X on the P*(0), the optimum control policies 
and the optimum trajectories. Consequently, the systems having initial 
values of X equal to 15, 10, and 5 were investigated. The P*(0) in each 
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case was obtained in a similar manner to that described for the system 
in which X(0) = 20. Optimal controls and optimal trajectories were 
obtained in each case. 

Upon investigation of the resulting optimal control policies, it 
was noted that they were remarkably similar in appearance' to the optimal 
control policy for the system having an initial X equal to 20. Fig. (4-5) 
shows the control policies for the four systems (X(0) = 20, 15, 10, 5). 
Note that the correspondence is very good for X(0) = 10 and 15 until the 
latter part of the common time interval at which point the three policies 
begin to separate. Of course, this correspondence only holds if the origin 
of the time axis for the systems having initial X equal to 10 and 15 are 
moved to the right as indicated. Notice, however, that the control policy 
for the system having an initial X equal to 5 does not seem to correspond 
at all. 

In spite of this latter result, the question arose as to whether the 
optimal control policy for the X(0) = 20 system, in conjunction with its 
trajectory, might provide the P*(0) for all initial values of X between 
20 and the lower limit of the trajectory of about 2.7. 

Consequently, systems having initial values of X of 17.5, 12.5, 7.5, 
2.5, and 1.0 were investigated for P*(0). Having obtained the P*(0) for 
these systems, the calculated P*(0) were ploted on a curve of optimal 
control versus optimal trajectory for the system having initial X equal 
to 20, Fig. (4-6). One immediately notices the close correspondence 
except for those system having initial X equal to 5, 2.5, and 1.0. Since 
the optimal trajectory for the X(0) = 20 system never got below a value 
of about 2.7, it is reasonable to assume that it would not be possible to 
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pondence except during the latter port of the common time interval and showing 
the non-correspondence for initial X equal to 5*0, 
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FIG. (4-6), Optimal Control plotted against Optimal Trajectory for the 
system having initial X equal to 20. Indicated points are the P*(o) obtained 
for systems having inital X equal to the values indicated. For T = 10 secs. 
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obtain a P*(0) for initial X which were less than this value, i.e., not 
on the curve. 

As a result of the above investigation, the possibility was suggested 
that if the system having initial X equal to 20 were allowed to run for 20 
seconds and a new optimal control policy and trajectory obtained for this 
system, that it might be that the control policy for the X(0) ^ 5 system 
would match and that this new U* vs, X* would provide the P*(0) for the 
values of X less than 2.7. 

The X(0) = 20 system was allowed to run for 20 seconds and a new 
optimum control and trajectory were obtained as well as a plot of U* vs. 
X*. When the P*(0) for the system investigated above were plotted on the 
new U* vs X* curve, Fig. (4-7), very close correspondence was noted in all 
cases including the ones for initial X of 5, 2.5 and 1.0. One must keep 
in mind that by extending the time interval to 20 seconds the problem 
became a new problem entirely. Even so, the results of this investiga- 
tion show that the U* vs. X* curve for X(0) = 20 over the extended time 
interval may be used to obtain an initial guess for the P*(0) for values 
of X(0) between 20 and zero. Modification of the initial guess in order 
to obtain the true P*(0) would be a much easier problem than the problem 
of obtaining P*(0) from scratch. 

The results, as stated above and as shown in the various figures, 
prove conclusively that it is possible to obtain numerical solutions to 
the problem of obtaining an optimal solution to the non-linear system by 
means of Pontryagin's Maximum Principle. 

It is to be pointed out that a very distressing feature of this 
method of solution is the excessive amount of time involved in obtaining 
a satisfactory solution to the problem compared v ith the problem time 
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FIG. (4-7). Optimal Control plotted against Optimal Trajectory for the 
system having initial X equal to 20. The time interval was extended in 
this case to 20 seconds. Indicated points are the P*(o) obtained for the 
systems having initial X equal to the values indicated. 
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of 10 or 20 seconds. (This would stggest some investigation into the 
possibility that more sophisticated programming of the CDC 1604 digital 
computer might result in a shorter solution time on the computer.) 

Finally, one additional fact should be borne in mind. Although the 
digital computer was utilized successfully to cbtain the optimal control 
function (specified at intervals of 0.1 seconds), building a physical 
component to duplicate this control is quite another matter. Some com- 
promise would probably be necessary. Perhaps a square pulse of U (as in 
the output of a zero order hold in a sample data system) corresponding in 
width or modulus to the peaks shown in Fig. (4-5) would suffice for engine- 
ering purposes. Another possibility would be to feed some percentage of 
the value of the trajectory back to the controller (negative feedback). 

In any case, each type of control would have to be evaluated and 
the cost function compared to the optimum one. It then becomes a 
question of deciding whether or not the results are sufficiently "opti- 
mum" for the engineering purpose in mind. Perhaps the only utility of an 
investigation such as this is in providing the "ideal” with which the 
engineer may judge the performance of the system he has designed. 

(Additional graphs of the optimal control policies and the optimal 
trajectories for the above mentioned initial X supplementing and support- 
ing the above data and conclusions may be found in Appendicies I and II.) 
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5. The cost function J 

The second system investigated was one in which the integral 



dt 



(5-1) 



was to be minimized with the same object equation as in Section 4. 



X = AX + B Sin X + CU 



(5-2) 



Proceeding in the same manner as in Section 4, introduce the vari- 



ables 



S 1 = X(t) 



■ t* 



dt , 



Substituting these variables in to Eqs. (5-1) and (5-2), and using 
the differential form yields the system equations 



X 1 = AX 1 + B Sin X x + CU 



x 2 = u 2 . 



By forming the final value functional 



S = x 2 (t) 



A X i = <¥l + C 2 X 2 



(5-3) 

(5-4) 

(5-5) 



the C vector is obtained. In order to minimize the final state of the 
adjoint equations, recall that P^(T) must equal to - C^. Thus the 
following boundary conditions for the adjoint space are obtained 



P 1 (T) = -C x = 0 

p 2 (t) = -c 2 = -1 

Forming the H- function as before 

H = P 1 (AX 1 + B Sin X 1 + CU) + 



(5-6) 



(5-7) 
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and performing the partial differentiations as defined by Eqs. (3-9), 



the following set of simultaneous equations is obtained, 

X x = AX^ + B Sin X 2 + CU 

X -u 2 

A 2 " U (5-8) 

P x = -(A + B Cos Xj^ 

P 2 = ° 

By differentiating H with respect to U and setting equal to zero, 
a relation between the control U, and the adjoint variables is obtained, 

CP, 



U = 



2P„ 



(5-9) 



Integrating the last of Eqs. (5-8), one is able to evaluate P£(t) 
as before. 



And since 



?2 = Constant. 



P 9 (T) = - 1 



then 

P 2 (t) = - 1 

Now substituting values for and C into Eqn. (5-9), the relation 
between P^ and U becomes 



U 



8 



(5-9a) 



Now Eqs. (5-8) may be arranged in a form suitable for solution 



X x = -.1X X + Sin X x + .03125 P x 


(5-10) 


Xg = P^ / 64 


(5-11) 


^ = (.lXj - Cos X 1 )P 1 


(5-12) 
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having the following boundary conditions^ 



X.^0) = 20 



P 1 (T) = 0 



(5-13) 



x 2 (o) = 0 



P 2 (T) = -1 



As in Section 4, two of the above equations, Eqs, (5-10) and (5-12), 
must be solved simultaneously after finding the P*(0) which minimizes the 
cost function, Eqn. (5-1). But, before stating the results of this in- 
vestigation, recall that the integral to be minimized is I dt 6 

Jo 

Even without Pontryagin's Maximum Principle, the solution is immediately 
obvious. U must equal zero if the cost function is to be minimized, since 
the control U, is independent of X in the object equation, Eqn. (5-1). 

Thus, one might hope for a reliable check on the answers obtained by the 
methods utilized in this paper in the solution to this system. 

Proceeding as in Section 4, values of P(0) between - 250 and + 50 
were investigated and the corresponding cost functions computed. The 
results verified the expected minimum of zero for P*(Q) = 0.0. Fig. (5-1) 
shows the minimum cost function for an initial P equal to zero which 
resulted in a zero control for the entire time interval; i.e. the un- 
controlled system. As was stated above, this is a reasonable and expect- 
ed solution. 

Utilizing the criteria J = \ I^dt provided a check on the solutions 

h 

obtained in synthesizing an optimal system by the maximum principle. 

The results obtained from the solution of Eqs c (5-10), (5-11) and (5-12) 
verified the predictions made through inspection of the H- function and 
its derivative, Eqs. (5-7) and (5-9). 

It appears from the above discussion, that this is not a very sensible 




as a cost function. A more meaningful use would be to 
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FIG. (5-1). Cost Function vs. Initial Value of P. X(0) - 20. 
j = f( u2)dt, showing a minimum cost function J - 0, at P(0) - 0 

apply the above cost function in a problem where the state of the tra 

jectory at time (T) = 10 secs, is fixed, say X(T) = 0.0. 

That is, to even have a reasonable problem we must have a fixed 

right end trajectory and also fix time, (T). 
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6. The cost function J 

The third system to be investigated was one in which the integral 



J = j 1 X 2 dt 

was to be minimized, again with the object equation 

X = AX + BSinT + CU 



(6-1) 



(6-2) 



Proceeding as before, it is seen that the system equations are 

(6-3) 



X x = AX X + B Sin X 1 + CU 



X 2 = ^ 



Forming the H-function 

H(P,X,U,t) = P 1 (AX 1 + B Sin X x + CU) + P^ 



(6-4) 



and taking suitable differentials, Eqs. (3-9), it is possible to obtain 
the adjoint equations. Four simultaneous equations are thus obtained 
which, when solved, provide the solution to the problem of obtaining 
the optimal system. 



X x = AXi + B Sin X 1 + CU 
X 2 = ^ 



P x = - P X (A + B Cos X x ) + 2?^ 

P 2 = ° 

From the final value functional, Eqn. (3-5), 

S =H C i X i = C l X l + C 2 X 2= £ X 1 

the C vector is obtained. As before, by setting P^T) equal to -C^ 
the terminal boundary conditions for the P^(t) are obtained. 



(6-5) 



dt 



(6-6) 
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(6-7) 



P 1 (T) = 0 



P 2 (T) = -1 



Since the last of Eqs. (6-5) may be integrated, it is seen that 



P 2 (t) = CONSTANT 



( 6-8 ) 



and, since 



P 2 (T) = -1 



it is seen that 



p 2 (t) = -i 



(6-9) 



Substituting this value for into the H- function, and then dif- 
ferentiating with respect to U and setting the resultant relation equal 
to zero in order to find the minimum, results in 



Since C is a constant equal to .25, the obvious conclusion must be that 



At first the last results seem quite disheartening as not relation 
between U and the adjoint space exist. However, this might have been 
anticipated since the H- function is linear in U and, as such, could have 
a minimum only at one of the boundaries. 

Reconsidering the Principle of the Maximum in view of the above 
results, one realizes that in order to minimize the H function, Eqn. 
(6-4), for any X and P, U must be as large as possible and of opposite 
sign to X (since the constant C is positive and it is desired to make 
that term negative). But, since most practical systems have an upper 




(6-10) 




( 6 - 11 ) 
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maximum on the amount of control available, U would then have to assume 
the largest value possible under such a constraint.* In mathematical 
notation, the above result would be expressed as 



u = ' x Mmx 



(6-12) 



As a result of the above developments, the only system equation 
needed to obtain the optimal control for this system is 



Eqn. (6-13) was solved numerically on the digital computer by means 
of the Runge-Kutta method for evaluation of Differential Equations. A 
computer solution was desired in order to verify the conclusions made 
on the basis of the above results and knowledge of the uncontrolled system, 
namely; (1) The cost function should decrease with increasing values of 
Control U, and (2), Once the system approached and was driven through 
zero, a chatter mode should be obtained for the remainder of the interval. 

The initial computer solutions immediately verified that the Cost 
Function did decrease as the absolute magnitude of the Control was in- 
creased. However, the chatter mode did not appear in the form expected. 
Fig. (6-1) shows that for a Control having a magnitude of 25, the Tra- 
jectory approaches zero in about 2,8 seconds and then begins to increase 
for a short while and then is driven back towards the zero point. This 
variation of the trajectory with time was not the chatter mode expected. 



*This result is similar to the well known Bang-Bang principle. 
Using Pontryagin's theory, this would be equivalent to minimizing Jdt 
used as the final value functional. 



X^ = AX^ + B Sin 



- S ^ X l |^| max 



(6-13) 
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FIG. (6-1). Trajectory, X(t), and Control, U(t) plotted as a func- 
tion of time, showing the chatter mode initially obtained from the 
solution of the Ja dt system. 

Note also that the control remained constant at -25 for the entire inter- 
val. The authors had expected that the negative control would drive the 
trajectory past the zero point, change sign as X(t) became negative ac- 
cording to Eqn. (6-12), and then drive X(t) positive. The cycle should 
repeat itself until the end of the interval. These results were not 
obtained as may be seen in Fig. (6-1). 

It was suggested that perhaps the sampling rate (0.1 sec) might be 
too coarse to show the chatter and that ’’noise" was interfering with the 
computations. Consequently, a sampling rate of .025 seconds was tried, 
but with no success. Results similar to those shown in Fig. (6-1) were 
obtained. 

At this point, the program being utilized in the solution of this 
system came under close scrutiny. It was discovered that the method 



34 



used to evaluate the differential equation (Runge-Kutta) utilized the 
last X(t) computed in the interval dt to set the sign of U(t) for the 
next evaluation of X(t) in the same interval. The program should have 
been using the sign of the X(t) computed during the n-lst time interval 
to control the sign of U(t) during the nth time interval. Because of 
the nature of the Runge-Kutta routine, four values of X(t) were being 
averaged in each time interval, some being positive values and some 
possibly negative, resulting in the erroneous results. 

The program was revised to correct the above deficiency and the 
results obtained showed the chatter mode and the alternating nature of 
the control as expected. See Fig. (6-2). 




FIG. (6-2). Sketch showing the chatter mode for X(t) and the 
alternating characteristics of U(t) for the system having AX^ dt 
as a cost function. 

The results obtained using the modified program confirmed the earlier 
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results that the cost function decreased with increasing magnitude of U. 




CONTROL |U| 

FIG. (6-3). A curve of Cost Function plotted against |uj showing 
the decrease in cost function (J) with the increase in |U| . 

Through investigation of this system employing various values of U, 
it was found that the magnitude of the control effected both the magni- 
tude and the frequency of the chatter. Additional graphical results 
may be found in Appendix IV. From an engineering point of view, it 
might be desirable to include a dead zone around zero in order to 
eliminate or reduce the chatter. 
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7 . Conclusions. 



In this investigation, the authors have inquired into and shown the 
feasibility of obtaining a numerical model for the optimal control of a 
non-linear system described by a first degree differential equation by 
means of the theory and procedures of Pontryagin's Maximum Principle as 
set forth by L. I. Rozonoer / 1 , 2, '/. The investigation was limited to 
the "free right end" problem and to the single dimensional case. 

The methods and procedures utilized may easily be extended to the 
n-dimensional case without any basic change in theory or procedure. 

Various cost functions have been investigated and some conclusions 
verified by numerical solution utilizing the CDC 1604 digital computer. 

Future investigation might well be in the area of second or third 
order systems and also in the area of computer programming techniques. 
The methods and programs utilized to obtain solutions were effective, 
but not very efficient. 

The theories utilized here do yield solutions. However, other 
methods (Dynamic Programming for instance) might well be more effective 
and useful to the investigator, depending upon his specific needs. 



37 



BIBLIOGRAPHY 



1. L. I. Rozonoer, "L. S. Pontryagin's Maximum Principle in the Theory 
of Optimum Systems -I", Automat ika i Telemakhanika, Vol 20, 1959, 

pp. 1320-1334 (English Translation in Automation and Remote Control) . 

2. L. I. Rozonoer, M L. S. Pontryagin's Maximum Principle in the Theory 
of Optimum Systems-II M , Automatika i Telemakhanika, Vol 20, 1959, 

pp. 1441-1458 (English Translation in Automation and Remote Control). 

3. L. I. Rozonoer, "L. S. Pontryagin's Maximum Principle in the Theory 
of Optimum Systems-III", Automatika i Telemakhanika, Vol 20, 1959, 
pp. 1561-1578 (English Translation in Automation and Remote Control). 

4. V. G. Boltyanskii, R. V. Gamkrelidge, and L. S. Pontryagin, "On the 
Theory of Optimum Processes", Doklady Akad. Nauk SSSR, Vol 110, No 1 
1956. (In Russian). 

5. R. Bellman, "Adaptive Control Processes: A Guided Tour, The RAND 

Corporation, 1961. 



38 



APPENDIX I 

GRAPHS FOR THE UNCONTROLLED SYSTEM 
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FIG. (i-l) Phase plane plot showing velocity as the ordinate against 
distance as the abscissa for the uncontrolled system having initial X equal 
to 20. This figure shows the graphical transformation of the abscissa from 
real time to distance. (Scale! Y-.05V^i ne » X-. 236 rad/line.) 
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FIG. (1-2) Uncontrolled System Trajectory (X) vs Time, for 11 positive 
values of X(o) showing the two stable equilibrium points, X = 2.85, and 842 
for the positive system. 
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Time 



(sec) 




FIG. (l-3) Uncontrolled System Trajectory (x) vs Time, for 11 negative 
values of X(0) showing the two stable equilibrium points, X =-2.85, and -8.52 
for the negative system. 
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APPENDIX II 

GRAPHS FOR THE COST FUNCTION J = 



j T (X 2 +U 2 )dt 
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GRAPHS FOR THE COST FUNCTION J = 



j 



(X" 

0 



U 2 ) 



dt 
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" = 1.57 

" = 0.218 

" = 0.024 

" = 2.77 

P*(0) = -128.0, 

(integer accuracy) 

J = 1156.01 and P*(0) 

= -127.428 8 J = 1099.52 

Same conditions as in 
II-8 above 

P*(0) = -93.82 

P*(0) = -41.96 

P*(0) = -14.7688 

Time interval = 20 sec. 
Time interval = 20 sec. 




local minima for the system having an initial X = 20, 
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FIG. (lI-2) Cost Function, J =jf(X*+ U 2 ’)dt, plotted as a function of 
initial value of the adjoint variable P, showing the true minimum and one 
local minimum for the system having an initial X = 17.5# 
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PIG. (II- 3 ) Cost Function, J =/(xS U z )dt, plotted as a function of 
initial value of the adjoint variable P, showing the true minimum and one 
local minimum for the system having an initial X _ 15. 
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FIG. (II- 4 ) Cost Function, J = /(X* + U 2 )dt, plotted as a function of 
initial value of the adjoint variable*?, showing the true minimum and two 
local minima for the system having an initial X = 7.5. 
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PIG. (II-6) 
initial value of 
local minima for 



Cost Function, J = Jj[x 2 + U z )dt, plotted as a function of 
the adjoint variable P, showing the true minimum and two 
the system having an initial X = 2.5. 




local minima for the system having an initial X = 1.0. 
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PIG. (lI-8) Optimal Trajectory (x) plotted against Time, for two 
values of P*(o) showing the variation in trajectory as a function of the 
accuracy to which P*(0) is determined; for the system having the Cost 
Function J = £(X*+ U z )dt and initial X = 20. 
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PIG, ( II— 10) Curves of Control (u) and Trajectory (x), plotted 
against Time, for the optimal system having the Cost Function J = ](x Z + U^)dt 
and initial X = 15. 0 
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Time (sec) 




PIG. (il-ll) Curves of Control (u) and Trajectory (x), plotted 
against Time, for the optimal system having the Cost Function J = J (X*+ U^)dt 
and initial X = 10. 
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PIG, (II-12) Curves of Control (u) amd Trajectory (x), plotted 
against Time, for the optimal system having the Cost Function J = j(X*+ U 2 )dt 
and initial X = 5. ' 
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FIG. (il- 13) Optimal Trajectory (x) plotted against Time, for the 
system were T final was extended to 20 seconds, having the Cost Function 
J = J(X 2 + U z ’)dt and initial X = 20.0. (Note when the time was extended 
from 10 to 20 seconds that the trajectory approached zero and was driven 
through the lower stable equilibrium point, X = 2.85 by the optimum controller J 
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FIG. (II-14) Optimal Control (u) plotted against Time, for the 
system were T final was extended to 20 seconds, having the Cost Function 
J =J 0 (X 2 + U l )dt and initial X = 20. 
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APPENDIX III 

GRAPHS FOR THE COST FUNCTION 




dt 
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APPENDIX IV 

GRAPHS FOR THE COST FUNCTION J 
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FIG. (iV-l) Curves of Control (u) and Trajectory (x), plotted against 
Time, for the optimal system having the Cost Function J = J (x 2 ) dt 
initial X = 20, and (Mote the absence of a chatter mode, as 

|U| wkK was not sufficient to drive the trajectory to zero in the alloted time.) 
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FIG. (lV-2) Optimal Trajectory (x) plotted against Time r for the 
system having the Cost Function J = J (X z )dt end initial X = 20. (Note the 
establishment of a chatter mode.) = 10. 
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FIG.(lV-3) Optimal Control (u) plotted against Time, for the system 
having the Cost Function J = j (Z^)dt and initial X = 20. (Note the establish- 
ment of a chatter mode.) !U| m ^ = 10» 
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PIG. (IV-4) Optimal Trajectory^/x) plotted against Time, for the 
system having the Cost Function J =X(X*)dt and initial X = 20. (Note the 
establishment of a chatter mode.) |U| M4x = 25. 
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PIG. (lV-6 ) Optimal Trajectory (x) plotted against Time, for the 
system having the Cost Function J = J,(x z )dt and initial X = 20. (Note the 
establishment of a chatter mode.) 1111^=50* 
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rra. (lV-7) Optimal Control (n) plotted against Tito, for the system 
having the Co.it Function J J t (:;''lt end initial X = 20. (Note the establish- 
raent of a chatter mode.) I^i***- >0. 
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APPENDIX Va 

SUBROUTINE RUNGE-KUTTA 
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V.a. The Runge-Kutta method of evaluating differential equations. 

There are many variants of the Runge-Kutta method* but the most 
widely used one is the following: given 

y' = f (x>y) 

y(x n ) = y n 

we compute in turn 

k l = h f( VV 

k. = h f(x + h/2, y + k /2) 

2 n n 1 

k = h f(x + h/2, y + k./2) 

3 n ■'n 2 

k 4 = h f(x n + h., y n + k 3 ) 

y n+l “ y n + 6 (k l + 2k 2 + 2k 3 + V • 

This process may be described in geometric terms. At the point 

(x n ,y n ) we compute the slope (k^/h) and using it we go one half step 
forward and examine the slope there. Using this new slope (k^/h) we 
again start at ( x n »y n )> 8° one half step forward, and again sample the 
slope. Using this latest slope (k^/h) we again start at ( x n »y. n ) but 
this time we go a full step forward where we examine the slope (k^/h) . 
The four slopes are averaged, using weights 1/6, 2/6, 2/6, 1/6, and 
using this average slope we make the final step from ( x n »y r ) to (x ^ , 

y ). If f(x,y) did not depend on y, then the averaging would produce 

5 

Simpsons' s formula. The method has an error term proportional to h , 

It is evident that the method throws away all old information and 
begins each complete step anew, and hence is hardly likely to be as 
efficient as methods which take advantage of old information. It is 
also evident that there is no check on whether the step size is too 
small or too large, though perhaps a study of the k^ might give a clue; 
this is not usually done. 
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The general spirit of the derivation is that the functions f(x,y) 
which are on the right-hand sides are all expanded in series in powers of 
h and the corresponding derivatives are equated to eliminate the lower 
powers of h. 

The method as used in a subroutine for the investigations of this 
thesis is given below in Fortran language. 



SUBROUTINE RKUTTA(N, T, X,DT) 

DIMENSION X(300), AK(4,300), XD0T(300), C(4) 

C(1)=0.0 

C (2)=0. 5 

C (3)=0. 5 

C (4)=1.0 

DO 4 1=1,4 

TC=T + C(I)*DT 

DO 2 J=1 , N 

2 XC(J)=X(J) + C(I)*AK(I-1), J) 

CALL DERIV (TC, XC, XDOT, N) 

DO 4 J=1,N 

4 AK(I, J)=DT*XDOT(J) 

DO 3 J=1,N 

3 X(J)+X(J) + (AK(1,J) +2.*AK(3, J) + AK(4,J))/6. 
RETURN 

END 
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APPENDIX Vb 

FORTRAN LANGUAGE PROGRAM TO OBTAIN TRAJECTORIES FOR THE 
UNCONTROLLED SYSTEM 
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no 



• • JOB GIBSON 9 JAN 1963 SPECIAL 
PROGRAM UNCON 

DIMENSION X ( 30 ) * Xl(900), 



10 



IU l«L VALUE, IT = rilNAL 

DT = TIME STEP, X(J) = ARRAY OF DEPENDENT VARIABLES 



It Y5 ( 900 ) , Y6 ( 900 ) . 

2Y11 (900), Y 1 2 ( 900 ) 
NUMPTS = 0 

READ 3, N, TC, TF, DT, 
N = NO* OF EQNS.t TO 
DT = TIME STEP, X 
F0RMAT(I10/(8F10*0)) 

T = TO 
NUMPTS = NUMPTS + 1 
XI (NUMPTS ) = T 
Y 1 (NUMPTS) = X( 1 ) 

Y2 ( NUMPT S ) = X ( 2 ) 

Y3 ( NUMPTS ) = X ( 3 ) 

Y4 ( NUMPTS ) = X ( 4 ) 

Y5 ( NUMPTS ) = X ( 5 ) 

Y6 ( NUMPT S ) - X ( 6 ) 

Y7 ( NUMPTS ) = X ( 7 ) 

Y8 ( NUMPTS ) = X ( 8 ) 

Y9 { NUMPTS ) = X ( 9 ) 

Y 1 0 { NUMPTS ) = X( 10) 

Y 11 (NUMPTS ) = X( 11 ) 

Y 1 2 ( NUMPTS ) = X ( 1 2 ) 

I F ( ,T F - T - DT) 10,20, 
CALL GRAPH (NUMPTS, XI 
CALL GRAPH (NUMPTS, XI 



RUN 

Y 1 ( 900 ) , Y2 ( 900 ) , Y3(900), Y4(900) 



Y7 ( 900 ) , Y8 ( 900 ) , Y9(900), Y10(900), 



( X( J),J=1 , N ) 
'INITIAL VALUE, 



TF * FINAL VALUE 



CALL GRAPH 
CALL GRAPH 
CALL GRAPH 
CALL GRAPH 
CALL GRAPH 
CALL GRAPH 
CALL GRAPH 
CALL GRAPH 
CALL GRAPH 
CALL GRAPH 



(NUMPTS, XI 
(NUMPTS, XI, 
(NUMPTS, XI , 
(NUMPTS, XI , 

( NUMPTS, XI , 
(NUMPTS, XI, 
(NUMPTS, XI , 
(NUMPTS, XI , 
(NUMPTS, XI . 
(NUMPTS, XI, 



Y 1 2 , 


8) 


YU, 


8) 


Y1C, 


8) 


Y 9 , 


8) 


Y 8 , 


8) 


Y 7 , 


8) 


Y 6 , 


8) 


Y 5 , 


8) 


Y 4 , 


8) 


Y 3 , 


8) 


Y2 , 


8) 


Yl, 


8) 



104 FORMAT (12H0 TRAJECTORY /(6E16.9)) 



PRINT 


104, 


( 


Yl 


(I ) , 


I = 


1, 


NUMPTS) 


PRINT 


104, 


( 


Y2 


(I), 


I = 


1, 


NUMPTS ) 


PRINT 


104, 


( 


Y3 


( I ), 


I = 


1, 


NUMPTS ) 


PRINT 


104, 


( 


Y4 


(I), 


I = 


1 , 


NUMPTS) 


PRINT 


104, 


( 


Y5 


(t). 


I = 


1, 


NUMPTS) 


PRINT 


104, 


( 


Y6 


(I), 


I = 


1 , 


NUMPTS ) 


PRINT 


104, 


( 


Y7 


(I ), 


I = 


1, 


NUMPTS ) 


PRINT 


104, 


( 


Y8 


(I), 


I = 


1, 


NUMPTS ) 


PRINT 


104, 


( 


Y9 


(I), 


I = 


1, 


NUMPTS) 


PRINT 


104, 


( 


Y10(I), 


I = 


1, 


NUMPTS ) 


PRINT 


104, 


( 


Y 11 ( I ) , 


I = 


1, 


NUMPTS ) 


PRINT 


1 04, 


( 


Y 1 2 ( I ) , 


I = 


1, 


NUMPTS) 


PRINT 


105, 


NUMPTS 









I 10) 



20 



DT) 



105 FORMAT (10H0NUMPTS = 

STOP 

CALL RKUTT A ( N , T, 

T = T ♦ DT 
GO TO 1 
END 

SUBROUTINE RKUTT A ( N , T,X ,DT ) 

DIMENSION X( 30 ) , AK (4 ,30) , XDOT ( 30 ) , XC ( 30 ) ,C(4) 
C(1 )=0.0 
C(2)=0.5 
C ( 3 ) =0* 5 
C ( 4 ) = 1 , 0 
DO 4 1=1,4 
TC=T+C ( I )*DT 
DO 2 J=1 , N 

XC( J) =X( J) ♦ C( I ) *A K ( 1-1, J) 

CALL DERIV(TC,XC,XDOT) 

DO 4 J= 1 , N 
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AK(I,J) = 
00 3 J = 1 
X( J ) =X { J)+ 
RETURN 
END 

SUBROUTINE 
DIMENSION 
XDOT (1) = 
XDOT (2) = 
XDOT (3) * 
XDOT (4) = 
XDOT (5) = 
XDOT (6) * 
XDOT (7) » 
XDOT (8) * 
XDOT (9) » 
XDOT ( 10 



DT «XDCT 1 J) 

• N 

( AK( 1 , J)+2.*AK (2, J)+2.*AK(3, J)+AK(4, J ) > /6 • 



DERI V ( T ,X »X OOT ) 
XDCT ( 30 ) y X ( 30) 



XDOT ( 
XDOT ( 
END 
END 

12 



11 ) 

12 



— C • 1 * X ( 1 

-0.1 * X ( 2 

-0.1 * X ( 3 

-0.1 * X ( 4 

-0.1 * X ( 5 

-0.1 * X ( 6 

-0.1 * X ( 7 

-C. 1 * X ( 8 

-0.1 * X ( 9 ) + SINF ( X< 9) ) 

) * -0.1 * X( 10 ) ♦ S I NF ( X( ion 
= -0. 1 • X( 1 1 ) + S I NF ( X( 1 1 > ) M 
> * -0.1 * X( 12 ) + S I NF ( X( 12 ) ) 



>+ SINF ( X( 1)) 
)+ SINF ( X( 2) ) 
)+ SINF ( X( 3)) 
>♦ SINF ( X( 4) ) 
)+ SINF ( X( 5) ) 
)+ SINF ( X( 6)) 
)+ SINF ( X( 7) ) 
)+ SINF ( X( 8)) 



.0 10.0 0.1 -.5 

-7.0 -7.5 -8.0 -8.5 

0 

02 X VS TIME UNCONTROLLED SYSTEM 
02 GIBSON AND ALMSTEDT 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 



-1.0 -2.0 -3. 

-10.0 -15.0 -20.0 

10 _ 08 

7 DECEMBER 1962 



-5.5 

2 
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APPENDIX Vc 

REPRESENTATIVE FORTRAN LANGUAGE PROGRAM TO OBTAIN COST 
FUNCTION AS A FUNCTION OF INITIAL P (ADJOINT VARIABLE) 
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non 



..JOB 



TIME 15 



300 



301 



302 



12 

21 

10 



11 

13 

101 

102 

20 



500 



GIBSON 9 JAN 1963 SPECIAL RUN MAX 
PROGRAM PJUX 
DIMENSION X ( 950 ) « Xl(950), Yl(950) 

READ 3, TO, TF, DT, XI, AL, AH, DA 

T0= INITIAL VALUE, TF= FINAL VALUE, DT = TIME STEP, XI* INITIAL X 
AH= HIGHEST INITIAL ADJOINT, AL* LOWEST INITIAL ADJOINT, 

D A= SPACING CF INITIAL ADJCINTS POINTS 
FORMAT (8F10.0) 

N = 900 

FORMS MATRIX OF INITIAL VALUES OF X 
DO 300 I * 1 ,N,3 
X ( I ) * XI 

FORMS MATRIX CF P-ZERO 
X ( 2 ) = AL 
DO 301 I * 5 ,N, 3 
X(I) = X( 1-3) + DA 

ZEROIZES INITIAL VALUES OF COST FUNCTION 
DO 302 I * 3, N, 3 
X(I) * 0.0 
T = TO 

NUMPTS * N/3 
L * 2 

DO 12 J = 1, NUMPTS 
X 1 ( J ) * X ( L ) 

L = L*3 

IF (TF -T) 1C, 20, 20 

DO ^1 J* 1 , NLMPTS 
Y 1 ( J ) = X C L ) 

L = L + 3 

DO 13 1=1, NUMPTS 

J = 2* ( I- 1 ) ♦ I 
XII) = X C J ) 

PRINT 101 

FORMAT ( 43H0P-ZER0 VS COST FUNCTION 
PRINT 102 ( X 1 ( I ) , Yl(I). X(I), 1 = 1, NUMPTS) 

FORMAT (F8.2, E16.9, E16.9) 

CALL GRAPH (NUMPTS, XI, Yl, 8) 

STOP 

CALL RKUTTA (N,T,X,DT) 

T = T + DT 
GO TO 21 
END 

SUBROUTINE RKUTT A( N, T ,X ,DT ) 

DIMENSION X ( 550 ) , AK(4,9$0), X00K950), XC(950), C(4) 

C( 1 )=0.0 
C ( 2 )=0.5 
C ( 3 ) =0. 5 
C ( 4 ) = 1 • 0 
DO 4 1=1 ,4 
TC = T +C ( I ) *DT 
DO 2 J* 1 , N 

XC ( J ) *X ( J ) + C(I)*AK(I— 1,J) 

CALL DERIV ( TC , XC, XDOT, N) 

DO 4 J* 1 , N 
AK ( I , J) * DT *XDOT ( J ) 

DO 3 J = 1 ,N 

X(J)=X(J) + (AK(1 , J)+2.*AK (2, J J+2.*AK(3, J)+AK(4, J ) )/6. 

RETURN 
END 

SUBROUTINE DER I V ( T ,X ,XDOT, N ) 

DIMENSION XDCT ( 950 ) , X(950) 

DO 500 K= 1 , N , 3 
L = K+ 1 
M=K+2 

XDOT ( K) = -0.1* X ( K ) + S IN F ( X ( K ) ) + .03125* X(L) 

XDO T ( L ) = 0.1* X ( L ) — X ( L ) * COSF(X(K))+ 2.0* X(K) 

XDOT ( M ) = X ( K ) * *2 + .015625* X(L)**2 
END 
END 



FINAL X) 



Tho following are the values of the variables read in for this program: 



TO = 0.0 
TF = 10.0 
DT = .1 
XI = 20.0 



AL = -250.0 
AH = 50.0 
DA = 1.0 
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APPENDIX Vd 

REPRESENTATIVE FORTRAN LANGUAGE PROGRAM USED TO COMPUTE 
THE OPTIMUM CONTROL POLICY AND TRAJECTORY HAVING THE 
P*(0) FOR THE SYSTEM. 
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■ JOB ALMSTEDT MK2 MOO 1 1/20/63 MAX TIME 5 MIN X(0) = 20.0 OPT. TF(20) 
PROGRAM MARK2 MOD 1 

7 DIMENSION X(30>, X 1 ( 10,201 ), Uf 1 0,201 ) , CFCN(IO), . 

7 1 XRA Y (201 ) t UNCLE (201 ) , TIME(201) 

X 1 = TRAJECTORIES, U= CONTROL EFFORT, CFCN= COST FUNCTION 
READ 101, N, TO, TF, DT, XI, AL 

N= TOTAL NUMBER OF EONS, TO = INITIAL TIME, TF = FINAL TIME 
XI = INITIAL VALUE OF X, AL = INITIAL P VALUE 
DO 1 I = 1 , N , 3 

1 X ( I ) = XI 

READ 102, (XU), 1=2, N, 3) 

DO 5 1= 3 , N, 3 

5 X(I)=0.0 
T * TO 
KK= 1 
NN = N/3 

30 I F ( TF - T) 10, 20, 20 

10 L = 3 

DO 1 1 1= 1 , NN 
CFCN ( I ) = X ( L ) 

11 L=L + 3 

T IME ( 1 ) = 0.0 
DO 15 I = 2, KK 

15 T I ME ( I ) = TIMEU-1M* DT 

DO 16 1= 1, NN 

PRINT 103, C FCN { I ) 

DO 17 K = 1, KK 
XRA Y ( K ) = X 1 ( I , K ) 

17 UNCLE (K) = U ( I , K) 

PRINT 1 04 , ( XRA Y (K), K = 1. KK ) 

PRINT 105, (UNCLE (K), K = 1,KK) 

NPTS = KK 

CALL GRAPH ( NPTS, TIME, XRAY , 8) 

CALL GRAPH ( NPTS, TIME, UNCLE, 8) 

16 CALL GRAPH ( NPTS, XRAY , UNCLE, 8) 

STOP 

20 L = 1 

DO 21 1= 1,NN 

X 1 ( I , KK ) = X ( L ) 

U ( I , KK ) = X(L+l)/8.0 

21 L= L+3 
KK= KK+ 1 

CALL RKUTT A ( N, T , X, DT ) 

T = T + DT 
GO TO 30 

101 FORMAT ( 110, 7F10.3) 

102 FORMAT ( 8F 1 0. 6) 

103 FORM AT ( 1 7 HOC OS T FUNCTION = E16.9)) 

104 FORMAT( 1 2H0T RA JECTORY /(6E16.9) 

105 FORMAT! 16H0C0NTR0L EFFORT /(6E16.9)) 

END 

SUBROUTINE RKUTT A ( N , T , X , DT ) 

DIMENSION X ( 60 1 ) , AK(4,601), XD0T(60l), XC(60l), C(4) 

C( 1 )=0.0 
C ( 2 ) =0. 5 
C ( 3 ) =0. 5 
C ( 4 ) =1 • 0 
DO 4 1=1,4 
TC=T+C ( I ) *DT 
DO 2 J= 1 , N 

2 XC(J) =X(J) + C(I)*AK(I-1,J) 

CALL DERI V (TC, XC, XDOT, N) 

DO 4 J=1 ,N 

4 AK ( I , J ) = DT *XDOT ( J ) 

DO 3 J = 1,N 

3 X( J)=X( JJMAKU , J)+2.*AM2, J)+2.*AK(3, J)+AK(4, J) )/6. 

RETURN 

END 

SUBROUTINE DERI V ( T, X, XDOT, N) 

DIMENSION XDOT (601 ) , X(601) 
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DO 500 K=1 ,N,3 
L =K+ 1 
M = K + 2 

XDOT(K)* -0.1* X ( K )♦ SINF(X ( K ) ) ♦ .03125* X(L) 
XDOT(L) * 0.1* X ( L ) - X( L ) * COSF(X(K)>+ 2.0* X(K) 
500 XDOT(M)* X ( K ) **2 ♦ .015625* X(L)**2 
END 
END 



90.0 
128.955976 



20.0 



0.1 



20.0 



X VS TIME 
ALMSTEDT 



P-ZERO - 
JAN 1963 



0.0 

10 



08 

20.0 NO 1 TF ( 20 ) 
MK2 ? 



10 08 
TF( 20 ) 



TRAJECTORY 
GIBSON AND 
0 

CONTROL FUNCTION U VS TIME NO 1 

GIBSON AND ALMSTEDT JAN 1963 MK2 

0 

10 08 

INITIAL CONTROL VS. TRAJECTORY OPTIMAL SYSTEM TF(20) 
GIBSON AND ALMSTEDT JAN 1960 MK2 

0 
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APPENDIX Ve 

FORTRAN LANGUAGE PROGRAM UTILIZED IN THE EVALUATION 
OF THE TRAJECTORIES AND CONTROL POLICIES FOR THE 
SYSTEM HAVING J(X 2 )dt AS THE COST FUNCTION. 



79 



‘..JOB GIBSON 22 JANUARY 1963 KK4MCD0 ALT 1 

PROGRAM MK4 MCDO ALT 1 

C THIS PROGRAM USES INTERGAL X SQUARE AS THE COST FUNCTION 

DIMENSION X ( 30 ) » XRAY( 1 5 ,500 ) , U(15), UNCLE ( 1 5 , 500 ) , TIME(500), 
1 Y1 ( 101 ) , Y2( 101 ) , V( 15) 

READ 101, NN, TO, TF, DT, XI, ( U ( I ) , 1=1 ,NN ) 

101 FORMAT (18/ 15F5.0) 

N = 2* NN 

C FORM X-MATRIX 

DO 1 1=1, NN 

1 X(I)=XI 
N 1 =NN + 1 

DO 2 1= N1 ,N 

2 X ( I ) = 0.0 
T = TO 
NPTS= 0 

27 I F ( TF - T) 10, 20, 20 

20 NPTS = NPTS + 1 

IF ( 1 - NPTS) 21,26,26 

21 CALL RKUTTA (N,T,X,DT, U, V) 

T = T + DT 

26 TIME(NPTS) = T 
DO 28 I = 1, NN 
IF ( X ( I ) ) 23, 24, 25 

23 UNCLE ( I , NPT S) = -U( I ) 

GO TO 22 

24 UNCLE ( I, NPTS) = 0.0 
GO TO 22 

25 UNCLE ( I , NPTS ) = U ( I ) 

22 XRAY ( I ,NP TS )= X ( I ) 

28 V ( I ) = UNCLEt I ,NPTS) 

GO TO 27 

10 DO 11 I = 1 , NN 
PRINT 111 , X ( 1 + N N ) 

111 FORMAT ( 1 7H0C0ST FUNCTION = E16.9) 

PRINT 112, (XRAY (I , J ) , J= 1 , NPTS ) 

112 FORMAT ( 12H0 TRAJECTORY / (6E16.9)) 

PRINT 113 (UNCLE ( I , J ) , J= 1,NPTS) 

113 FORMAT ( 1 6H0CGNTR0L EFFORT /(6E16.9)) 

DO 12 K= 1 , NPTS 

Y1 (K) = XR AY ( I , K ) 

12 Y2(K) = UNCLEt I,K) 

CALL GRAPH ( NPTS , TIME , Yl , 8 ) 

11 CALL GRAPH ( NP TS , T I ME , Y2 , 8 ) 

PRINT 114, NPTS 

114 FORMAT ( 8H0NPTS = 15) 

STOP 



RKUTTA (N,T,X,DT, U, V) 

X ( 30 ) , AK ( 4 ,3 0 ) , XDOT (30), XC(30), C(4),V(15), U(30) 



END 

SUBROUTINE 
DIMENSION 
C( 1 )=0.0 
C ( 2 ) =0 . 5 
C ( 3 ) =0. 5 
C ( 4 ) = 1 . 0 
DO 4 I = 1 ,4 
TC = T + C ( I ) *DT 
DO 2 J =1 , N 

XC ( J ) = X ( J ) + C ( I ) * AK ( I - 1 , J ) 

CALL DERI V ( XC, TC,* XDOT, N, U, V) 

DO 4 J = 1,N 
A K ( I , J ) = DT *XDGT ( J ) 

DO 3 J = 1 ,N 

X(J)=X(J)+(AK( 1»J)+2.*AK(2,J)+2.*AK(3,J)+AK(4,J) )/6. 

RETURN 

END 

SUBROUTINE DERIV ( X , T , XDOT , N, U , V) 

DIMENSION XCOT ( 30 ) , X(30), U ( 1 5 ) , V ( 1 5 ) 

NN = N/2 
Nl = NN + 1 
DO 1.1= 1 , NN 

XDOT ( I ) = -.1*X(I) + S INF ( X ( I ) ) +.25* V(I) 

DO 2 I = N1 ,N 

XDOT (l)= X ( I-NN ) **2 

END 

END 



0.0 



10 

10 . 0.1 



20.0 -5.0-10.-15. -20. -25. -30. -35. -40. -45. -50. 
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